Using FEMLAB for Gravitational problems: numerical simulations for all. 
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We discuss the possibility to solve Modern Numerical Relativity problems using finite element 
methods (FEM). Adopting a "user friendly" software for handling totally general systems of nonlin- 
ear partial differential equations, FEMLAB, we model and numerically solve in a short time a Gowdy 
vacuum spacetime, representing an inhomogeneous cosmology. Results agree perfectly with existing 
simulations in the recent literature based not of FEMs but on finite differences methods. Possible 
applications for non relativistic Astrophysics, General Relativity, elementary particle physics and 
more general theories of gravitation like EMDA and branes are discussed. 
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I. INTRODUCTION 



General Relativity (GR) deals with dynamical defor- 
mations of space and time, massively using, as it happens 
for continuum mechanics, tensor calculus. Mathemati- 
cally, General Relativity is described by a very large num- 
ber of coupled strongly nonlinear PDEs. The extraordi- 
nary developments of Numerical Relativity of the last ten 
years, due to rapid developments of large scale computa- 
tional resources, have shown that some tools developed 
in the past by mechanical engineers, and in particular 
Finite Elements techniques, can be fruitfully adopted in 
this field too. The spontaneous formation of event hori- 
zons, Cauchy horizons and infinite curvature singulari- 
ties, still poses serious problems both at theoretical and 
numerical level. For this reason Numerical Relativity is 
today a branch of theoretical physics which is growing up 
in complexity, requiring knowledge of both high level GR 
theoretical background as well as purely numerical one. 
Numerics in particular still represents a terrible obstacle 
for those scientists which want to enter such a stimu- 
lating new area. In this article we show that FEMLAB 
p , a user friendly software developed for solving general 
systems of nonlinear partial differential equations from 
zero (ODEs) up to three dimensions, with or without 
time dependence, in general as well as in weak form, can 
be usefully applied to gravitational physics problems. In 
Q, FEMLAB was used for studying the scattering on 
a sonic black hole, applying standard procedures devel- 
oped for analyzing wave fields around Kerr black hole, 
i.e. a 3 + 1 constrained evolution scheme and excision in 
horizon penetrating coordinates 3j. Although the prob- 
lem was linear, the presence of an event horizon required 
an extremely accurate use of FEMLAB due to necessity 
of very high resolution and precision in order to numeri- 
cally keep all the perturbations confined inside the black 
hole once they have crossed its event horizon and avoid 
constraint violations. Here instead, we attack an exact 
nonlinear problem, i.e. we study a 1 + 1 set of Einstein's 
vacuum field equations describing a Gowdy inhomoge- 



nous cosmology, showing how to model in less then ten 
minutes the problem with FEMLAB. Other ten minutes 
of running on a standard laptop pc are required in or- 
der to be ready to compare the results of the simulation 
with existing literature based on standard and advanced 
finite differences methods. It appears clear that FEM- 
LAB is able to solve perfectly the problem capturing the 
fine structure (the spikes) characteristic of this class of 
universes. This application suggests that FEMLAB can 
become for sure a necessary tool for those scientists and 
mathematicians working in the field of full GR analy- 
sis as well as in partial differential equations present in 
non relativistic astrophysics, elementary particle physics, 
statistical mechanics as well as extended theories of grav- 
itation (Brans-Dicke, branes, EMDA, etc.). 



II. FIELD EQUATIONS 

In this section we introduce the mathematical frame- 
work for our FEMLAB simulations. We follow reference 
y (and references therein) for derivation of the field 
equations and comparison of numerical solutions. The 
Gowdy spacetime on T 3 x R has the form 
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where P, Q and A are functions of t and x. The T 3 spa- 
tial topology is imposed by having < x,y,z < 2ir and 
being P, Q and A periodic functions of x. We introduce 
the coordinate r = — In t such that the singularity is ap- 
proached as r — > oo. In terms of this coordinate, vacuum 
Einstein field equations split into first order evolution 
equations in divergence form (comma denotes ordinary 
derivation) 

P,r=R, 
Q ,r = S , 
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which trivially determine A once P and Q are known. 
For this reason we will not be interested in this second 
set of equations which does not influence the evolution 
set once a proper initial data satisfying the constraints is 
chosen. In particular, a meaningful initial data at r = 
is P = 0, R = vo cos x : Q = cos x, S = (where vq = 5 in 
our simulations). Periodic boundary conditions are im- 
posed at x = and x — 2tt, i.e. P(t, 0) = P(t, 27r), t > 
and similarly for and 5, while "cusps" are avoided 
by imposing zero flux Neumann boundary conditions at 
both ends, i.e. Q yX = P iX = only (in fact R^ x and 
5 ja ; = do not appear in the field equations © for the 
evolution). This is a delicate point and must be clarified, 
because periodicity would require to equate gradients at 
both ends without specifying their value (as we did in our 
case of zero flux instead). Initial data however is sym- 
metric with respect to x = 7r and field equations, naively, 
"do not prefer" right or left directions (change x <-> — x 
in equations (0) and @), consequently evolution must 
be symmetric on both sides. Matching solutions at both 
ends, we can avoid cusps only if gradient is zero there (no 
flux). Clearly this naive statement can be formalized rig- 
orously by using the evolution equations (J2J cast in null 
coordinates and characteristic variables 3 . Due to these 
symmetries, the simple zero flux conditions generate a 
smooth matching on both ends, so in this sense the initial 
request of periodicity for (P, Q,R,S) has become useless. 
We will continue to require such a condition only to show 
how FEMLAB can handle periodic boundary conditions. 
We can now start our finite element modelling, pointing 
out that an excellent introduction to the involved math- 
ematical theory of Finite Elements can be found in the 
first chapters of ref.0 
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FIG. 1: Choosing the model properties. 
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FIG. 2: Drawing the spatial domain. 



III. BUILDING THE MODEL 

In figure we start selecting the dimensionality of 
the problem (ID), the type of problem (time dependent 
in general form), the variables (P, Q, R, S) and the ele- 
ment (Lagrange quintic). We adopted the highest non- 
linear element in order to get the best results, although 
a quadratic Lagrange element leads practically the same 
results (with smaller computational time). In figure (0 
we draw the domain, i.e. a line starting at x = and 
ending at x = 2tt. In figures ©-© we insert the field 
equations in divergence form while in figure © we enter 
the initial data. In these FEMLAB modeling pictures 
as well as in the forthcoming text and figures, time t has 
always to be intended as scaled time r of equations J2J). 
In figures J7J) and (0 instead we chose zero flux Neumann 



boundary conditions at both ends of the domain, while 
in figure one enters the periodic boundary conditions 
(which as described before, are automatically implied by 
zero flux ones). We show one figure only because the en- 
tire procedure would take several slides due to the fact 
that one has to make the value of the various field co- 
incide at both ends (vertices). The procedure is very 
simple, although takes few minutes and is well explained 
in the online manual with several examples so we refer 
to it for this part of the model. In figures (|10 )) -(|T2 J) 
we select the numerics of the problem. A direct solver 
(UMFPACK) is chosen, time integration starts at t = 
and ends at t = 10, saving on hard disk the simulation 
results at every 0.1 time step. Other options are chosen 
in order to have the most precise settings for the sim- 
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FIG. 3: Entering field equations (first part). 
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FIG. 6: Entering initial data. 
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FIG. 4: Entering field equations (second part). 
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FIG. 5: Entering field equations (third part). 



FIG. 7: Selecting zero flux conditions (first part). 
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FIG. 8: Selecting zero flux conditions (second part). 
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FIG. 9: Selecting periodic boundary conditions (see text). 
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FIG. 10: Selecting numeric properties (first part). 



ulation. Time stepping is automatically chosen by the 
solver which refines it when necessary (although one can 
always set it up manually). Concerning numerical ac- 
curacy, we have selected relative and absolute tolerance 
thresholds of 10 -7 . This choice will be sufficient to get 
high quality results. In figure ((13)) we select the mesh- 
ing of the domain, which is uniform and has a spacing of 
Ax = 0.005. While FEMLAB can handle adaptive mesh 
refinement, we prefer to use such a very fine constantly 
spaced meshing here. 



IV. RESULTS 
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FIG. 11: Selecting numeric properties (second part). 
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FIG. 12: Selecting numeric properties (third part). 



This visualization clearly shows the generation and sub- 
sequent development of these Gowdy spikes. The initial 
data for Q, i.e. cosx is not evident in figure ([T%|) due to 
the very large change of scale, which has passed from 1 
to around 100 in ten time units. 



As anticipated before, simulation takes approximately 
ten minutes to run on a Pentium 4 running at 1.7GHz. In 
figure {H} we show P at t = 10, in the interval < x < 
while figure ([T5|) is a plot of Q with same conditions. Fig- 
ure (|16|) shows a detail (a spike) of figure ((14)). These 
three figures are exactly coincident with figures 1,2, and 
3 of reference 4], i.e. FEMLAB produced the same re- 
sults of finite differences simulations using finite elements 
instead. In figures (fTT)) and (|TH|) we show the time evolu- 
tion in the entire spatial domain of P and Q respectively. 



V. DISCUSSION 

The simulation presented in this article and the com- 
parison with existing literature suggests that FEMLAB 
can become for sure a necessary tool for those scientists 
and mathematicians working in the field of GR analysis, 
as in the past years Maple with its tensor packages did. 
Due to its possibility to handle general PDEs, problems 
in astrophysics (equilibrium configurations, MHD, ...) 
and in more general theories of gravitation like EMDA 
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FIG. 16: Detail of the plot for P at t — 10 (read r) around 
x 0.36. 

FIG. 13: Selecting meshing properties. 




FIG. 15: Plot of Q at t = 10 (read r) for < x < tt. FIG - 18: Time evolution of Q in the interval < x < 2tt. 
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and branes can be solved. Being FEMLAB totally inte- 
grated in MATLAB, and having MAPLE a direct export 
to it, the process of symbolic derivation of field equation, 
numerical solution and export for deeper analysis ( for 
studying the solution in the frequency domain in order 
to extract quasi- normal modes, for example) and com- 
parison with experimental data results straightforward. 
Moreover, being FEMLAB designed for 64bit platforms, 
large amounts of memory (typically a ten of Gb) can be 
addressed for simulations in 2D and 3D, although for the 
latter case, having in mind General Relativity problems 
like strongly distorted black holes and neutron stars, a 
parallelized version of the software would certainly be 
welcome. Even for those working in Numerical Relativ- 
ity, the use of FEMLAB allows them to test an idea in a 



simplified version in few hours, and if the obtained results 
appear to be promising, write down then an optimized 
code for their large scale supercomputers. 
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